﻿using System;
using System.Linq;

namespace Algorithm_Comparison
{
    class LHD
    {
        public readonly int NumberOfPoints; // number of points of this LHD
        public readonly int Dimension; // dimension of this LHD
        public readonly double[,] DistanceMatrix; // matrix containing all the distances
        public readonly double[] PointMinDistance; // distance to the closest point
        private readonly Point[] _points; // array containing all the points of this LHD

        public double MinDistance; // distance between the 2 closest points
        public double LastCalculatedPhiP; // So we don't have to recalculate PhiP all the time

        public LHD (int numberOfPoints, int dimension, Point[] points)
        {
            NumberOfPoints = numberOfPoints;
            Dimension = dimension;
            DistanceMatrix = new double[NumberOfPoints, NumberOfPoints];
            PointMinDistance = new double[NumberOfPoints];

            //Initialize points
            _points = points;
            for(int pointIndex = 0; pointIndex < NumberOfPoints; pointIndex++)
            {
                _points[pointIndex] = new Point(Dimension);
            }
        }

        public LHD(int numberOfPoints, int dimension) : this(numberOfPoints, dimension, new Point[numberOfPoints]) { }

        public LHD(LHD original) : this(original.NumberOfPoints, original.Dimension)
        {
            for(int pointIndex = 0; pointIndex < NumberOfPoints; pointIndex++)
                _points[pointIndex] = new Point(original._points[pointIndex]); // copy the points from lhd:original
            DistanceMatrix = (double[,]) original.DistanceMatrix.Clone(); // copy the distance matrix
            PointMinDistance = (double[]) original.PointMinDistance.Clone(); // copy the point min distances
            MinDistance = original.MinDistance;
            LastCalculatedPhiP = original.LastCalculatedPhiP;
        }
        /// <summary>
        /// Randomizes the coordinates of the points of this LHD
        /// Range for coordinates is [0, NumberOfPoints-1]
        /// </summary>
        public void SetRandomMatrixILS()
        {
            for (int columnIndex = 0; columnIndex < Dimension; columnIndex++)
            {
                int[] randomPerm = GenerateRandomPermutation(NumberOfPoints);

                for (int pointIndex = 0; pointIndex < NumberOfPoints; pointIndex++)
                {
                    _points[pointIndex].Coordinates[columnIndex] = randomPerm[pointIndex];
                }
            }
        }
        /// <summary>
        /// Randomizes the coordinates of the points of this LHD
        /// Range for coordinates is [0, (NumberOfPoints-1)/NumberOfPoints]
        /// </summary>
        public void SetRandomMatrixESE()
        {
            for (int columnIndex = 0; columnIndex < Dimension; columnIndex++)
            {
                int[] randomPerm = GenerateRandomPermutation(NumberOfPoints);

                for (int pointIndex = 0; pointIndex < NumberOfPoints; pointIndex++)
                {
                    _points[pointIndex].Coordinates[columnIndex] = ((double)randomPerm[pointIndex])/NumberOfPoints;
                }
            }
        }
        private static readonly Random Generator = new Random();
        /// <summary>
        /// Generates a random permutation
        /// </summary>
        /// <returns> A random permutation of 0..numberOfPoints-1 </returns>
        public static int[] GenerateRandomPermutation(int numberOfPoints)
        {
            var result = new int[numberOfPoints];
            for (int i = 0; i < numberOfPoints; i++)
                result[i] = i;

            for (int i = 0; i < numberOfPoints; i++) // randomize the sequence
            {

                int next = Generator.Next(i, numberOfPoints);
                int buffer = result[i];
                result[i] = result[next];
                result[next] = buffer;
            }
            return result;
        }

        /// <summary>
        /// Recalculates the entire distance matrix
        /// </summary>
        /// <param name="euclidian">Whether or not to use euclidian distances</param>
        public void CalculateDistanceMatrix(bool euclidian=true)
        {
            for (int pointIndex = 0; pointIndex < NumberOfPoints; pointIndex++)
            {
                for (int pointIndex2 = pointIndex + 1; pointIndex2 < NumberOfPoints; pointIndex2++)
                {
                    DistanceMatrix[pointIndex2, pointIndex] = _points[pointIndex].DistanceTo(_points[pointIndex2], euclidian);
                    DistanceMatrix[pointIndex, pointIndex2] = DistanceMatrix[pointIndex2, pointIndex];
                }
            }
        }

        /// <summary>
        /// Recalculates only the corresponding row and column
        /// </summary>
        /// <param name="pointIndex">The column and row to recalculate</param>
        /// <param name="euclidian"></param>
        public void CalculateDistanceMatrix(int pointIndex, bool euclidian=true)
        {
            for (int pointIndex2 = 0; pointIndex2 < NumberOfPoints; pointIndex2++)
            {
                if (pointIndex2 != pointIndex)
                {
                    DistanceMatrix[pointIndex2, pointIndex] = _points[pointIndex].DistanceTo(_points[pointIndex2], euclidian);
                    DistanceMatrix[pointIndex, pointIndex2] = DistanceMatrix[pointIndex2, pointIndex];
                }
            }
        }

        /// <summary>
        /// Calculates the distance to the closest point for each point (distance matrix needs to be correct)
        /// Also calculates the distance between the 2 closest points (CalculateMinDistance())
        /// </summary>
        public void CalculatePointMinDistance()
        {
            
            for (int pointIndex = 0; pointIndex < NumberOfPoints; pointIndex++)
            {
                double smallestDistance = double.MaxValue;
                for (int pointIndex2 = 0; pointIndex2 < NumberOfPoints; pointIndex2++)
                {
                    if (pointIndex2 != pointIndex)
                    {
                        double interDistance = DistanceMatrix[pointIndex, pointIndex2];
                        if (smallestDistance > interDistance)
                            smallestDistance = interDistance;
                    }
                }
                PointMinDistance[pointIndex] = smallestDistance;
            }
            CalculateMinDistance();
        }

        public void CalculateMinDistance()
        {
            MinDistance = Int32.MaxValue;
            foreach (double pointMinDist in PointMinDistance)
                if (pointMinDist < MinDistance)
                    MinDistance = pointMinDist;

        }

        /// <summary>
        /// Returns the number of points of which the distance to the closest point is equal to distance
        /// </summary>
        /// <param name="distance"></param>
        /// <returns></returns>
        public int MinDistanceOccurences(double distance)
        {
            return PointMinDistance.Count(dist => dist == distance);
        }

        /// <summary>
        /// Returns the number of points for which the distance to the closest point is equal the MinDistance
        /// </summary>
        /// <returns></returns>
        public int MinDistanceOccurences()
        {
            return PointMinDistance.Count(dist => dist == MinDistance);
        }

        /// <summary>
        /// True if the distance to the closest point of the point with index==pointIndex is equal to the smallest distance (MinDistance)
        /// </summary>
        /// <param name="pointIndex"></param>
        /// <returns></returns>
        public bool IsCriticalPoint(int pointIndex)
        {
            return PointMinDistance[pointIndex] == MinDistance;
        }

        /// <summary>
        /// Returns a new latin hyper cube of which the coordinate value of point1 and point2 is swapped in column columnIndex
        /// Also updates the distance matrix and the pointMinDistances and the MinDistance of the new latin hyper cube
        /// </summary>
        /// <param name="point1"></param>
        /// <param name="point2"></param>
        /// <param name="columnIndex"></param>
        /// <param name="euclidian"></param>
        /// <returns></returns>
        public LHD LocalMove(int point1, int point2, int columnIndex, bool euclidian=true)
        {
            LHD lhdResult = new LHD(this);
            double coordPoint1 = lhdResult._points[point1].Coordinates[columnIndex];
            double coordPoint2 = lhdResult._points[point2].Coordinates[columnIndex];

            lhdResult._points[point1].Coordinates[columnIndex] = coordPoint2;
            lhdResult._points[point2].Coordinates[columnIndex] = coordPoint1;

            lhdResult.CalculateDistanceMatrix(point1, euclidian);
            lhdResult.CalculateDistanceMatrix(point2, euclidian);

            lhdResult.CalculatePointMinDistance();
            return lhdResult;
        }

        /// <summary>
        /// Returns whether or not this latin hyper cube is better than "other" based on the Maximin criterion
        /// </summary>
        /// <param name="otherLHD"></param>
        /// <returns></returns>
        public bool IsBetterMaximinThan(LHD otherLHD)
        {
            if (otherLHD == null)
                return true;
            return MinDistance > otherLHD.MinDistance ||
                   (MinDistance == otherLHD.MinDistance && MinDistanceOccurences() < otherLHD.MinDistanceOccurences());
        }

        /// <summary>
        /// Performs a cyclic order exchange in column columnIndex from pointIndex1 to pointIndex2
        /// Aso updates the distance matrix, the point min distances and the mindistance
        /// </summary>
        /// <param name="pointIndex1"></param>
        /// <param name="pointIndex2"></param>
        /// <param name="columnIndex"></param>
        public void CyclicOrderExchange(int pointIndex1, int pointIndex2, int columnIndex)
        {

            double lastPointCoordinate = _points[pointIndex2].Coordinates[columnIndex];

            for (int pointIndex = pointIndex2; pointIndex > pointIndex1; pointIndex--)
                _points[pointIndex].Coordinates[columnIndex] = _points[pointIndex-1].Coordinates[columnIndex];

            _points[pointIndex1].Coordinates[columnIndex] = lastPointCoordinate;

            for (int pointIndex = pointIndex1; pointIndex <= pointIndex2; pointIndex++)
                CalculateDistanceMatrix(pointIndex);
            CalculateMinDistance();
        }

        /// <summary>
        /// calculates the PhiP value
        /// </summary>
        /// <param name="pValue"></param>
        /// <returns></returns>
        public double CalculatePhiP(int pValue) // old method
        {
            double result = 0d;
            foreach (double d in PointMinDistance.Distinct().OrderBy(d => d))
                result += MinDistanceOccurences(d) / Math.Pow(d,pValue);

            result = Math.Pow(result, 1d / pValue);
            LastCalculatedPhiP = result;
            return result;
        }
        //Alternative method which didn't work for some reason
        /**public double CalculatePhiP(int pValue)
        {
            double result = 0d;
            for (int i = 0; i < NumberOfPoints; i++)
                for (int j = i + 1; j < NumberOfPoints; j++)
                    result += 1d/Math.Pow(DistanceMatrix[i, j], pValue);

            result = Math.Pow(result, 1d / pValue);
            LastCalculatedPhiP = result;
            return result;
        }/**/

        /// <summary>
        /// Calculates the function value which will be used in ESEMaximin
        /// </summary>
        /// <returns></returns>
        public double CalculateMaximinESE()
        {
            return MinDistanceOccurences(MinDistance) / MinDistance;
        }

        ///<summary>
        /// Converts the coordinates of all the points from ESE format to ILS format
        /// basically multiplies all the coordinates by NumberOfPoints
        /// </summary>
        public void ConvertPointsFromEseToIls()
        {
            foreach (Point p in _points)
            {
                for (int column = 0; column < Dimension; column++)
                    p.Coordinates[column] *= NumberOfPoints;
            }
        }
        /// <summary>
        /// Converts the coordinates of all the points from ILS format to ESE format
        /// Basically divides al lthe coordinates by NumberOfPoints
        /// </summary>
        public void ConvertPointsFromILSToEse()
        {
            foreach (Point p in _points)
            {
                for (int column = 0; column < Dimension; column++)
                    p.Coordinates[column] /= NumberOfPoints;
            }
        }

    }
}
